
% Monte Carlo simulations 

function [v1_mean, v2_mean, alpha_reshape, alpha_hat_reshape, beta_M_reshape, beta_v1_reshape, beta_v2_reshape,S_reshape, SM, ER_reshape, ERM, Data]=Betas_2D(h0, H, G, Gamma, T, dt, r, N, v1_bar, kappa1, c1, v2_bar, kappa2, c2, mu_D, Strike, lambda, p, K, rho)

N_sq=N^2;
Nf=4; %Number of factors+1. When Nf=5 and additional factor V is included. 


M=round(T/dt); % number of random variables for similations

Data_S=zeros(K,N_sq);
Data_M=zeros(K,1);
Data_v1=zeros(K,1);
Data_v2=zeros(K,1);
Data_V=zeros(K,1);
X=ones(K,Nf); 
alpha=zeros(N_sq,1);
alpha_hat=zeros(N_sq,1);
beta_M=zeros(N_sq,1);
beta_v1=zeros(N_sq,1);
beta_v2=zeros(N_sq,1);
beta_v3=zeros(N_sq,1);
alpha_hat_fm=zeros(N_sq,1);
beta_M_fm=zeros(N_sq,1);
beta_v1_fm=zeros(N_sq,1);
beta_v2_fm=zeros(N_sq,1);
beta_v3_fm=zeros(N_sq,1);
q=zeros(N_sq,1);

%dw=randn([M,K]).*sqrt(dt);
%dw_isc=randn([M,K]).*sqrt(dt);


tic
xi_0=0;
F=0;
F0=0;
FM=0;
EM=0;  ES=0;
v0=abs(v1_bar)+0.01;


for k=1:1:K % index of a simulation
    rng('shuffle');
    Sigma1=0;
    Sigma2=0;
    V=0;
    v1=v1_bar;
    v2=v2_bar;
    rng shuffle; % creates a different seed each time
    dw=randn([M,1]).*sqrt(dt);
    rng shuffle; % creates a different seed each time
    dw_tilde=rho.*dw+sqrt(1-rho^2).*randn([M,1]).*sqrt(dt);
    for m=1:1:M % index of the time period
    % Process for state variables
        v1=max(0.0001,v1+kappa1*(v1_bar-v1)*dt+c1*sqrt(v1)*dw(m));% systematic volatility
        v2=max(0.0001,v2+kappa2*(v2_bar-v2)*dt+c2*sqrt(v2)*dw_tilde(m));% idiosyncratic volatility
        V=V+sqrt(v1)*dw(m);
        Sigma1=Sigma1+v1*dt;
        Sigma2=Sigma2+v2*dt;
    end 
    
    % Real option part of the aggregate consumption
    C_RO=exp(lambda*mu_D*T-0.5*lambda*Sigma1.*H.^2-0.5*lambda*(1-lambda).*G.^2.*Sigma2+lambda*V.*H)*Strike^(1-lambda);
    
    % Aggregate consumption
    C_T=p*N_sq*exp(mu_D*T-0.5*h0^2*Sigma1+h0*V)+(1-p)*sum(C_RO);
    % SPD at date T
    xi=C_T.^(-Gamma);
    F=F+xi.*exp(lambda*mu_D*T-0.5*lambda*Sigma1.*H.^2-0.5*lambda*(1-lambda).*G.^2.*Sigma2+lambda*V.*H)./K;
    F0=F0+xi*exp(mu_D*T-0.5*h0^2*Sigma1+h0*V)/K;
    FM=FM+xi.*C_T./K;
    % SPD at date 0
    xi_0=xi_0+xi*exp(r*T)/K;
    
    S_T=p*exp(mu_D*T-0.5*h0^2*Sigma1+h0*V)+(1-p).*Strike^(1-lambda).*exp(lambda*mu_D*T-0.5*lambda*Sigma1.*H.^2-0.5*lambda*(1-lambda).*G.^2.*Sigma2+lambda*V.*H);
    
    EM=EM+C_T./K; 
    ES=ES+S_T./K;%stock expected value
    
    Data_S(k,:)=S_T;
    Data_M(k,1)=C_T;
    Data_v1(k,1)=Sigma1;
    Data_v2(k,1)=Sigma2; 
    Data_V(k,1)=V; 
end

F=F./xi_0;
F0=F0/xi_0;
S=p*F0+(1-p)*Strike^(1-lambda).*F; % Stock price at date 0
SM=FM./xi_0; % Market price at date 0

% Calculating stock and market returns
for k=1:1:K
    Data_S(k,:)=Data_S(k,:)./S-1;
end
Data_M(:,1)=Data_M(:,1)./SM-1;

ER=(ES-S)./S;
ERM=(EM-SM)./SM;


% ICAPM alpha without factor mimicking

Data=[Data_S, Data_M, Data_v1, Data_v2];
X(:,2)=Data_M;
X(:,3)=Data_v1;
X(:,4)=Data_v2;
X(:,5)=Data_V; % in the case one additional factor is added
%Omega=(X(:,1:Nf)'*X(:,1:Nf))\X(:,1:Nf)';
Omega=inv(X(:,1:Nf)'*X(:,1:Nf))*X(:,1:Nf)';

for i=1:1:N_sq
    Y=Data(:,i)-r;
    vec=Omega*Y;
    alpha_hat(i)=vec(1);
    beta_M(i)=vec(2);
    beta_v1(i)=vec(3);
    beta_v2(i)=vec(4);
end

% CAPM alpha

Omega=inv(X(:,1:2)'*X(:,1:2))*X(:,1:2)';
e=ones(N_sq,1);
for i=1:1:N_sq
    Y=Data(:,i)-r;
    vec=Omega*Y;
    alpha(i)=vec(1);
    e(i)=i;
end
%-------------------------------------------------------------------------

% Factor mimicking portfolio and alpha

% Factor mimicking portfolio for systematic volatility



Matrix=[X(:,1) X(:,2) X(:,3)];
Omega=(Matrix'*Matrix)\Matrix';
q=(Omega*Data_S)';


% Sorting by slope coefficients
vec_sort=sortrows([q(:,3)  e]);
ind=round(N_sq/5);
ind1_vec=vec_sort(1:ind,2)';
ind2_vec=vec_sort(ind+1:2*ind,2)';
ind3_vec=vec_sort(2*ind+1:3*ind,2)';
ind4_vec=vec_sort(3*ind+1:4*ind,2)';
ind5_vec=vec_sort(4*ind+1:5*ind,2)';


VIX1=sum(Data_S(:, ind1_vec),2)/ind;
VIX2=sum(Data_S(:, ind2_vec),2)/ind;
VIX3=sum(Data_S(:, ind3_vec),2)/ind;
VIX4=sum(Data_S(:, ind4_vec),2)/ind;
VIX5=sum(Data_S(:, ind5_vec),2)/ind;

Matrix=[X(:,1) VIX1 VIX2 VIX3 VIX4 VIX5];

Omega=(Matrix'*Matrix)\Matrix';
vec1=Omega*Data_v1;
Data_v1_fm=Matrix*vec1-vec1(1); 
v1_mean=mean(Data_v1_fm);


% Factor mimicking portfolio for idiosyncratic volatility


Matrix=[X(:,1) X(:,2) X(:,4)];
Omega=(Matrix'*Matrix)\Matrix';
q=(Omega*Data_S)';


% We regress on Sigma_1 rather than h_i^2*Sigma_1, and then after
% estimating coefficient in front of Sigma_1, we do a correction 
% q(:,3)=q(:,3)./(G.^2)'; 

% Sorting by slope coefficients
vec_sort=sortrows([q(:,3)  e]);
ind=round(N_sq/5);
ind1_vec=vec_sort(1:ind,2)';
ind2_vec=vec_sort(ind+1:2*ind,2)';
ind3_vec=vec_sort(2*ind+1:3*ind,2)';
ind4_vec=vec_sort(3*ind+1:4*ind,2)';
ind5_vec=vec_sort(4*ind+1:5*ind,2)';

VIX1=sum(Data_S(:, ind1_vec),2)/ind;
VIX2=sum(Data_S(:, ind2_vec),2)/ind;
VIX3=sum(Data_S(:, ind3_vec),2)/ind;
VIX4=sum(Data_S(:, ind4_vec),2)/ind;
VIX5=sum(Data_S(:, ind5_vec),2)/ind;

Matrix=[X(:,1) VIX1 VIX2 VIX3 VIX4 VIX5];

Omega=(Matrix'*Matrix)\Matrix';
vec2=Omega*Data_v2;
Data_v2_fm=Matrix*vec2-vec2(1); 
v2_mean=mean(Data_v2_fm);

%==========================================================================
%--------------------------------------------------------------------------

X(:,2)=Data_M;
X(:,3)=Data_v1_fm;
X(:,4)=Data_v2_fm;
Omega=(X(:,1:Nf)'*X(:,1:Nf))\X(:,1:Nf)';

for i=1:1:N_sq
    Y=Data(:,i)-r;
    vec=Omega*Y;
    alpha_hat_fm(i)=vec(1);
    beta_M_fm(i)=vec(2);
    beta_v1_fm(i)=vec(3);
    beta_v2_fm(i)=vec(4);
end

alpha_hat_reshape=reshape(alpha_hat_fm, N, N);
beta_M_reshape=reshape(beta_M_fm, N, N);
alpha_reshape=reshape(alpha, N, N);
ER_reshape=reshape(ER, N, N);
S_reshape=reshape(S, N, N);

%--------------------------------------------------------------------------

% Systematic volatility beta
 X(:,2)=Data_M;
 X(:,3)=Data_v1_fm;
 Omega=(X(:,1:3)'*X(:,1:3))\X(:,1:3)';
 
 for i=1:1:N_sq
     Y=Data(:,i)-r;
     vec=Omega*Y;
     alpha_hat_fm(i)=vec(1);
     beta_v1_fm(i)=vec(3);
 end

beta_v1_reshape=reshape(beta_v1_fm, N, N);

% IVol beta

X(:,3)=Data_v2_fm;
Omega=(X(:,1:3)'*X(:,1:3))\X(:,1:3)';

for i=1:1:N_sq
    Y=Data(:,i)-r;
    vec=Omega*Y;
    alpha_hat_fm(i)=vec(1);
    beta_v2_fm(i)=vec(3);
end

beta_v2_reshape=reshape(beta_v2_fm, N, N);

toc
end
